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QUALITY  CONTROL  OF  METEOROLOGICAL  OBSERVATIONS 


Toke  Jayachandran  and  Richard  Franke 

Int roduct  ion :  A  major  input  component  to  weather  prediction 
systems  is  the  observational  data,  typically  geopotential  heights 
and  winds  from  weather  observing  stations  around  the  world.  The 
Navy  Operational  Global  Atmospheric  Prediction  System  (NOGAPS) 
produces  weather  forecasts  using  a  spectral  forecast  model. 
Preceding  the  forecast  is  a  multivariate  objective  analysis, 
followed  by  a  nonlinear  normal  mode  initialization  scheme.  The 
multivariate  objective  analysis  scheme  uses  the  weather 
observations  to  produce  analyzed  fields  with  minimum  statistical 
error.  As  is  to  be  expected,  some  of  the  observations  will  turn 
out  to  be  erroneous  for  various  reasons  such  as  the 
mis-calibrat ion  of  the  measurement  devices,  failure  to  perform 
prescribed  data  corrections,  communications  errors  and  simple 
recording  errors .  We  present  in  this  paper  a  quality  control 
methodology  to  screen  the  data  to  identify  those  observations  that 
appear  to  be  erroneous  in  a  statistical  sense.  It  is  then  up  to  a 
meteorological  analyst  or  a  decision  making  algorithm  to  examine 
the  tagged  observations  and  decide  on  a  course  of  action  -  delete 
the  miscreant  observations,  modify  the  data  if  appropriate  or 
include  the  observations,  as  is,  for  NOGAPS  analysis.  The 
procedures  can  also  be  used  to  monitor  the  reliability  of  the 
observing  stations  and  to  alert  them  of  any  mispract ices . 

Data  Quality  Tests:  The  methodology  proposed  here  is  an 
adaptation  of  a  statistical  test  for  outliers  (Dixon,  1950  [1]). 
The  first  step  in  the  application  of  this  test  is  to  partition  the 
mc.nitoring  stations  inlo  small  "contiguous  groups".  The  recorded 
observations  within  any  single  group  are  expected  to  be  "similar" 
because  the  weather  patterns  at  these  monitoring  sites  are 
comparable  due  to  the  proximity  of  the  stations  to  each  other. 
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One  approach  towards  identifying  such  a  group  of  stations  is  the 
following.  Select  the  latitude  and  longitude  for  a  target 
location  on  the  surface  of  the  earth  (this  could  be  the 

coordinates  of  a  reporting  station  of  interest)  and  choose  an 

approximately  square  region  by  specifying  a  range  for  the 
latitude,  say  5°  above  and  below  the  target,  and  within  a  similar 
distance  in  the  longitudinal  direction.  The  observational 
increments  (differences  between  the  reported  measurements  and  the 
forecast  values)  for  all  of  the  stations  within  the  designated 
area  are  then  subjected  to  a  sequence  of  outlier  tests. 

Assume  that  the  number  of  reporting  stations  within  the  target 
area  is  N  and  the  data  from  the  most  recent  )c  days  (say  30  days) 

is  to  be  examined  for  quality.  Let  Xu  ,  Xi2  i  •  •  •  >  be 

the  observation  increments  for  the'  ith  station,  i  =  1,  2,  .  .  .  , 

N;  let  Ml  ,  M2  ,  .  .  .  ,  M^  be  the  means  of  the  increments 

(averaged  over  the  k  days  for  each  of  the  N  stations)  and  Si,  S2  , 

.  .  .  ,  Sfj  the  standard  deviations  (over  the  k  days)  for  the  N 

stations . 

The  first  test  we  propose  is  to  compare  the  daily  observational 
increments  Xij  ,  X2j  ,  ,  Xnj  for  the  jth  day,  and  check  for 

"outliers".  Let  Yi  <  Y2  <  .  .  .  <  Yn  be  the  data  rearranged  in 

increasing  order.  Compute 

rii  =  (Yn  -  Yn-i)  /  (Yn  -  Y2)  ; 

i f  rii  exceeds  the  tabulated  value  in  Table  A  (extracted  from 
Dixon, 1950) ,  corresponding  to  a  chosen  statistical  significance 
level  li  (say  .05),  conclude  that  the  observed  datum  for  the 
station  corresponding  to  Yn  is  an  outlier.  This  implies  that  the 
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observation  increment  is  too  large  relative  to  the  increments 
for  the  other  stations  in  the  region.  That  is,  the  observed 
geopotential  height  appears  to  be  too  high  (large)  when  compared 
to  the  forecast .  Now  compute 


r'li  =  (Y2  -Yi)  /  (Yn-1  -  Yi) 

and  conclude  that  Yi  is  an  outlier  if  r’n  exceeds  the  tabulated 
value  in  Table  A;  if  this  happens,  the  implication  is  that  the 
observed  geopotential  height  is  too  low  (small)  compared  to  the 
forecast .  These  tests  can  be  performed  for  each  of  the  days  for 
which  data  is  available.  The  meteorological  analyst  or  the 
decision  making  algorithm  will  then  have  to  choose  a  course  of 
action,  based  on  some  established  decision  rule,  for  each  of  the 
observations  failing  the  test.  The  two  tests  were  applied  to 
actual  data  from  three  different  areas  of  the  globe  viz.,  US, 
Europe,  and  China.  In  each  area  we  identified  a  location  by 
specifying  its  latitude  and  longitude  and  then  selected  all  the 
stations  within  4°  of  this  location.  The  results  of  the  daily 
tests  are  presented  in  Tables  la,b,c.  The  I.D.  numbers  of  the 
stations  included  for  the  test  are  listed  in  the  first  column. 
The  entries  (there  are  two  entries  in  each  of  the  16  columns)  in 
the  table  represent  the  number  of  days  in  October  '90  on  which  the 
stations  failed  the  outlier  test  at  the  low  end  (rn'  test)  and  the 
test  at  the  high  end  (r^  test)  using  the  data  (geopotential  height 

differences)  for  16  pressure  levels  (1000,  850,  700,  500,  400, 
300,  250,  200,  150,  100,  70,  50,  30,  20,  15,  10  mb) .  A  dot  in  any 
position  implies  no  failures;  an  asterisk  indicates  insufficient 
data  to  perform  the  tests.  Such  a  table  can  be  used  in  making 
quality  judgments  on  the  performance  of  the  monitoring  stations. 

Another  test  that  can  also  be  used  to  monitor  the  performance 
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quality  of  weather  stations  is  to  perform  the  two  outlier  tests  on 
the  averages  ,  M2  ,  .  .  .  ^  .  The  data  for  the  tests  are 

these  averages  rearranged  in  increasing  order.  Tables  2a,b,c 
present  the  same  type  of  information  as  Tables  la,b,c  except  that 
this  time  the  entries  are  a  1  or  a  dot;  a  1  signifies  that  the 
average  for  a  station  failed  the  outlier  test  while  a  dot  would 
signify  that  the  station  passed  the  test. 

A  third  outlier  test,  that  can  be  viewed  as  a  test  for 
consistency,  is  to  check  for  an  outlier  in  the  station  standard 
deviations  Si  ,  S2  ,  .  .  .  ,  .  If  Yi  <  y2  <  .  .  .  <  Yn  are 
these  standard  deviations  in  ordered  form,  compute 

rio  =  (Yn  -  Yn-i)  /  (Yn  -  Yi) 

and  conclude  that  the  standard  deviation  corresponding  to  Yn  is 
statistically  too  large  relative  to  the  others,  if  rio  exceeds  the 
critical  value  in  Table  B  (extracted  from  Dixon,  1950)  at  a 
significance  level  ii  (say  .05).  Tables  3a,b,c  are  similar  to 
Tables  1  and  2  using  the  station  standard  deviations,  with  the 
difference  that  only  the  high  end  outlier  test  is  used. 

We  have  written  a  computer  program,  using  Fortran,  that  will 
select  (1)  the  reporting  stations  within  a  specified 
(approximately  square)  regior*  around  a  specified  location  and  (2) 
the  appropriate  data  for  these  stations  from  the  Naval 
Oceanographic  and  Atmospheric  Research  Laboratory  (NOARL)  files 
[2],  and  automatically  performs  all  of  the  tests  described  above. 
Tables  1,  2  and  3,  described  abovt,  are  the  actual  outputs  of  this 
program.  The  program  also  generates  additional  information  such 
as  the  number  of  times  the  daily  outlier  tests  were  performed 
(Table  4)  and  the  number  of  days  for  which  data  was  available 
(Table  5),  during  the  period  specified  for  the  test.  The  numbers 
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in  Table  4  will  sometimes  be  smaller  than  the  corresponding 
entries  in  Table  5;  this  happens  whenever  the  number  of  data 
points  available  for  the  test  is  smaller  than  the  minimum  sample 
size  of  4.  Table  6  is  an  example  of  another  type  of  output  of  the 
program;  for  a  selected  station  the  table  provides  the  daily 
test  history,  over  the  16  pressure  levels.  A  "PASS"  indicates 
that  the  station  passed  all  the  daily  tests;  a  "LFAIL"  signifies 
that  the  observational  increment  is  too  low  as  compared  to  the 
numbers  for  the  other  stations  in  the  region,  on  at  least  one  day 
of  the  testing  period.  A  "FAILH"  denotes  that  the  increment  for 
at  least  one  day  was  deemed  to  be  high.  Similar  tables  displaying 
the  results  of  the  means  and  standard  deviations  tests  are  also 
produced.  In  the  "TOTAL"  column  a  "PASS"  means  that  the  station 
passed  all  the  tests  on  all  the  days  and  a  "FAIL"  implies  the 
converse . 
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STATION  /  Lvl  1  2  3 


4  5 


6 


7  8  9  10  11  12  13  14  15  16 


72456 

72357 

72349 

72340 

72260 

72247 

72239 

72257 


• 

2 

1 

1 

1 

1 

3 

2 

2 

1 

1 

2 

* 

. 

2 

5 

2 

1 

1 

2 

• 

1 

1 

1 

2 

2 

1 

1 

* 

* 

• 

1 

1 

1 

* 

* 

• 

1 

1 

1 

* 

* 

• 

2 

1 

1 

2 

1 

1 

1 

2 

* 

* 

1 

1 

i 

i 

1 

• 

1 

, 

1 

1 

1 

3 

5 

* 

1 

• 
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* 

* 

* 

t 

* 

* 

* 

* 

* 

* 

♦ 

* 

* 

* 

* 

* 

* 

* 

♦ 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

t 

* 

* 

♦ 

* 

t 

* 

* 

* 

* 

* 

* 

Table  la:  Failures  in  daily  outlier  test  (.  =0,  *  =  no  test) 
U.S.  stations  near  (35°, -95°) 

10/01/90  to  10/31/90  at  OOZ 


STATION  /  Lvl  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16 


26702 

4 

1  . 

1 

.  1 

.  3  . 

1  . 

•  0  •  » 

•  •  • 

.  1  . 

1 

.  1 

1 

2 

1 

26406 

1 

.  1 

2  . 

1.1, 

.  1 

.  .  3  . 

.  .  1 

, 

1 

12425 

.  1 

• 

.  1 

«  •  • 

.  1  .  . 

1  .  . 

1 

.  2 

i 

2 

i 

1 

12374 

•  ■ 

.  1  .  : 

1  . 

•  •  «  • 

1  .  1 

.  1  . 

1 

*  , 

• 

2 

2 

1 

12330 

.  1 

1 

«  • 

.  .  .  1  . 

•  •  •  • 

*  •  • 

.  1  . 

• 

1 

1 

* 

* 

* 

* 

12120 

, 

»  • 

.  .  1 

■  • 

•  •  •  * 

•  •  • 

•  •  • 

• 

* 

* 

* 

10437 

i  . 

1 

«  • 

•  «  «  1  • 

•  •  •  • 

«  •  • 

*  * 

* 

* 

4c 

* 

* 

* 

* 

10338 

«  • 

•  •  « 

•  •  •  • 

0  0  0 

•  *  • 

0  0 

0 

0 

, 

0 

0 

0 

10046 

i 

.  i 

2  . 

•  •  •  « 

0  0  0 

•  •  • 

*  * 

* 

* 

* 

* 

* 

* 

* 

* 

10035 

•  • 

•  t  «  • 

0  0  0 

•  •  • 

0 

• 

• 

1 

• 

0 

9393 

, 

«  • 

.  .  1 

•  •  «  • 

0  0  0 

•  •  • 

,  , 

0 

• 

1 

, 

• 

0 

9184 

* 

«  • 

•  «  « 

•  •  •  « 

.  1  . 

.  .  1 

0 

» 

* 

• 

0 

6181 

•  • 

•  •  • 

•  •  •  • 

•  *  • 

•  •  • 

0 

• 

0 

2591 

•  • 

.  .  1 

2 

.1.1 

•  •  • 

•  • 

• 

• 

* 

* 

* 

* 

2527 

•  • 

.  .  1 

1 

•  •  •  • 

•  «  • 

•  •  • 

•  • 

, 

♦ 

* 

* 

* 

Table  lb:  Failures  in  daily  outlier  test  (,  =0,  *  =  no  test) 
European  stations  near  (55  ,15°) 

10/01/90  to  10/31/90  at  OOZ 


STATION  /  Lvl  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16 


57679 

1 

1 

.  1  . 

2 

* 

* 

57516 

* 

* 

1 

1 

1 

1  . 

1 

2 

1 

1 

1 

1 

* 

* 

57494 

1 

1 

2 

1 

1 

1 

i 

1 

* 

* 

57461 

* 

* 

* 

1 

1 

* 

* 

5744  7 

3 

* 

* 

* 

* 

57328 

* 

* 

* 

♦ 

1 

1 

* 

* 

* 

*  t 

* 

* 

* 

* 

* 

* 

♦ 

* 

* 

* 

* 

57178 

, 

* 

* 

i 

2 

.  1  1 

1 

1 

2 

3  . 

1 

1 

* 

* 

57127 

2 

1 

1 

1 

1  . 

1 

1 

* 

* 

57083 

1 

1 

1 

1 

* 

* 

57036 

4 

1 

i 

i 

1 

1 

1 

* 

* 

57245 

* 

* 

* 

* 

,  * 

♦ 

* 

* 

* 

* 

* 

* 

* 

* 

* 

57633 

* 

* 

* 

* 

* 

1 

.  * 

* 

* 

* 

* 

* 
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*  * 

* 

* 

t 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

Table  Ic:  Failures  in  daily  outlier  test  (.  =0,  *  =  no  test) 
Chinese  stations  near  (32°, 110°) 

10/01/90  to  10/31/90  at  12Z 
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STATION  /  Lvl  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16 


1  .  1 


. 1.1... 

t***************tt************** 
ifttttt**********  **************** 


Table  2a:  Failures  in  the  means  test  (.  =0,  *  =  no  test) 
U.S.  stations  near  (35°, -95°) 

10/01/90  to  10/31/90  at  OOZ 


STATION  /  Lvl  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16 


26702  . 

26406  . 

12425  . 

12374  . *  * 

12330  . *  *  *  * 

12120  . ■ . ****** 

10437  .  ********** 

10338  . 

10046  .  ********** 


10035 

9393 

9184 

6181 


2591  . **** 

2527  . **** 

Table  2b:  Failures  in  the  means  test  (.  =0,  *  =  no  test) 


European  stations  near  (55°, 15°) 
10/01/90  to  10/31/90  at  OOZ 


72456 

72357 

72349 

72340 

72260 

72247 

72239 

72257 


STATION  /  Lvl  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16 


57679  . ♦  * 

57516  **.l . 1 . 1.1.1.. 

57494  . 

57461  .,** . ** 

57447  .  .  *  *  . 

57328  **** . 1.  .  ,***^************** 

57178  .  .  *  * . *  * 

57127  . 

57083  . *  * 

57036  . *  * 

57245  .  ********** 

57633  ******************************** 


Table  2c:  Failures  in  the  means  test  (.  =0,  *  =  no  test) 
Chinese  stations  near  (32  ,110  ) 

10/01/90  to  10/31/90  at  OOZ 
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STATION  /  Lvl  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16 

72456  . 

72357  . 

72349  . 

72340  . 

72260  . 

72247  .  1  ...  1  1  . 

72239  **************** 

72257  *************:*:♦* 

Table  3a:  Failures  in  the  standard  deviations  test  (.  =0,  *  =  no  test) 

U.S.  stations  near  (35°, -95°) 

10/01/90  to  10/31/90  at  OOZ 


STATION  /  Lvl  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16 

26702  . 

26406  . 

12425  .  1  .  .  . 

12374  * 

12330  .1  . *  * 

12120  . ; . *  *  * 

10437  ***** 

10338  . 

10046  ***** 

10035  . 

9393  . 

9184  . 

6181  . 

2591  *  * 

2527  *  * 

Table  3b:  Failures  in  the  standard  deviations  test  (.  =0,  *  =  no  test) 

European  stations  near  (55°, 15°) 

10/01/90  to  10/31/90  at  OOZ 


STATION  /  Lvl  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16 

57679  . * 

57516  »  . 

57494  . 

57461  .  * . * 

57447  .  *  . 

57328  ** . ********* 

57178  .  * . * 

57127  . 

57083  . * 

57036  . * 

57245  ** . ***** 

57633  **************** 

Table  3c:  Failures  in  the  standard  deviations  test  (.  =0,  *  =  no  test) 

Chinese  stations  near  (32  ,110°) 

10/01/90  to  10/31/90  at  OOZ 
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STATION 

/  Lvl  1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

72456 

31 

30 

31 

31 

31 

31 

31 

31 

31 

31 

30 

29 

27 

24 

23 

0 

72357 

31 

30 

31 

31 

31 

30 

31 

31 

31 

30 

30 

30 

29 

26 

22 

0 

72349 

31 

30 

31 

31 

30 

31 

30 

30 

30 

30 

30 

29 

28 

28 

22 

0 

72340 

31 

29 

31 

31 

30 

30 

30 

30 

30 

30 

29 

28 

28 

26 

22 

0 

72260 

30 

29 

30 

30 

30 

30 

30 

29 

29 

29 

29 

29 

28 

26 

23 

0 

72247 

31 

28 

31 

31 

31 

31 

31 

31 

31 

31 

29 

29 

29 

27 

23 

0 

72239 

1 

1 

1 

1 

1 

1 

1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

72257 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Table  4: 

Number  of 

days  daily  outlier 

test 

was 

performed 

U.S. 

stations 

near 

(35 

o 

1  “ 

■95°) 

10/01/90 

to 

10/31/90  at 

OOZ 

STATION  / 

Lvl  1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

72456 

31 

30 

31 

31 

31 

31 

31 

31 

31 

31 

30 

29 

28 

25 

24 

11 

72357 

31 

30 

31 

31 

31 

30 

31 

31 

31 

30 

30 

30 

30 

28 

25 

4 

72349 

31 

30 

31 

31 

30 

31 

30 

30 

30 

30 

30 

29 

28 

28 

23 

6 

72340 

31 

29 

31 

31 

30 

30 

30 

30 

30 

30 

29 

28 

28 

27 

24 

2 

72260 

30 

29 

30 

30 

30 

30 

30 

29 

29 

29 

29 

29 

29 

28 

27 

3 

72247 

31 

28 

31 

31 

31 

31 

31 

31 

31 

31 

29 

29 

29 

27 

23 

6 

72239 

1 

1 

1 

1 

1 

1 

1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

72257 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Table  5:  Number  of  daily  reports 
U.S.  stations  near  (35*^, -95°) 
10/01/90  to  10/31/90  at  OOZ 


STATION  72456;  LAT.LON  =  39,  -95 


LEVEL 

#DAYS 

DAILY 

TESTS 

MEANS  TESTS 

STD  TESTS 

TOTAL 

1 

31 

PASS 

PASS 

PASS 

PASS 

PASS 

PASS 

2 

30 

PASS 

PASS 

PASS 

PASS 

PASS 

PASS 

3 

31 

PASS 

PASS 

PASS 

PASS 

PASS 

PASS 

4 

31 

PASS 

*FAILH 

PASS 

PASS 

PASS 

♦FAIL* 

5 

31 

PASS 

PASS 

PASS 

PASS 

PASS 

PASS 

6 

31 

LFAIL* 

PASS 

PASS 

PASS 

PASS 

♦FAIL* 

7 

31 

LFAIL* 

PASS 

PASS 

PASS 

PASS 

♦FAIL* 

8 

31 

LFAIL* 

*FAILH 

PASS 

PASS 

PASS 

♦FAIL* 

9 

31 

LFAIL* 

PASS 

PASS 

PASS 

PASS 

♦FAIL* 

10 

31 

PASS 

PASS 

PASS 

PASS 

PASS 

PASS 

11 

30 

LFAIL* 

PASS 

PASS 

PASS 

PASS 

♦FAIL* 

12 

29 

PASS 

PASS 

PASS 

PASS 

PASS 

PASS 

13 

28 

LFAIL* 

PASS 

PASS 

PASS 

PASS 

♦FAIL* 

14 

25 

PASS 

*FAILH 

PASS 

PASS 

PASS 

♦FAIL* 

15 

24 

LFAIL* 

*FAILH 

PASS 

PASS 

PASS 

♦FAIL* 

16 

11 

NOTEST 

NOTEST 

PASS 

PASS 

PASS 

PASS 

Table 

6 :  Summary 

report 

for  station 

72456 

10/01/90  to  10/31/90  at  OOZ 
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on  no 


Appendix:  Program  listing 

C  THIS  PROGRAM  IS  FOR  CHECKING  THE  STANDARD  FORECAST  DEVIATION 

C  FILES  FOR  OUTLIERS  IN:  (1)  DAILY  READINGS,  (2)  MEANS,  AND 

C  (3)  STANDARD  DEVIATIONS  OF  STATIONS  WITHIN  A  TARGET  DISTANCE 

C  OF  A  GIVEN  LATITUDE  AND  LONGITUDE. 

C 

C  WRITTEN  BY  RICHARD  FRANKE 

C  DEPARTMENT  OF  MATHEMATICS 

C  NAVAL  POSTGRADUATE  SCHOOL 

C  MONTEREY,  CA  93943 

C  0083P@CC.NPS.NAVY.MIL 

C  408/646-2758 

C 

C  DATE  OF  BEGINNING:  6/21/91 

C  DATES  OF  CHANGES:  6/21/91 

C  MODIFIED  TO  USE  SYMBOLS  FOR  OUTPUT:  6/21/91 

C  MINOR  PRETTYING  UP:  6/24/91 

C  MODIFIED  TO  NOTE  NUMBER  OF  TESTS  AND  SUMMARY  BY  STATION:  6/27/911 

C 

PARAMETER  ( MX=152 , NS=15 ) 

REAL  FE{MX,NS,16) ,FER(MX,NS),XMN{NS),STD(NS),FED(NS) 

I NTEGER  YRMD ( MX , NS ) , NSTAT ( NS ) , KNT ( NS ) , NSD ( NS ) , 

1  NDLF(NS,16) ,NDHF{NS,16) ,NSDF{NS,16),NMNLF(NS,16) ,NMNHF(NS,16) , 

2  NSS(NS) ,YMD(MX,NS) ,NDG(NS,16),NDDT(NS,16) ,LATS(NS),LONS(NS) 
CHARACTER*15  lUNIT.OUNIT 

CHARACTER*6  FLAGS (-1:1), FLAGSL ( -1 : 1 ) , FLAGSH (-1:1) 

CHARACTER*!  SY(-1:36) 

CHARACTER* 38  SYY 
EQUIVALENCE  (SV,SYY) 

D AT A  SIT/ ’ * . 1 2 34 5 6 7 8 9 ABCDE FGH I J KLMNOFQRSTU VWXY Z & ’ / 

DATA  FLAGS/’ NOTEST’ , ’  PASS  ’,’*FAIL*’/, 

1  FLAGSL/’ NOTEST’ , ’  PASS  ’,’LFAIL*’/, 

2  FLAGSH/’NOTEST’ , ’  PASS  ’,’*FA1LH’/ 

C 

GET  DATA  INPUT  FILE  NAME  AND  OUTPUT  FILE  NAME 
WRITE(*,9) 

9  FORMAT (’  NAME  THE  INPUT  FILE’) 

READ(*,10)IUNIT 
10  FORMAT! A15) 

OPEN(4,BLANK=’ZERO’ ,FILE=IUNIT,MODE=’READ’ ,STATUS=’OLD’ ) 

WRITE (*,21) 

21  FORMAT!’  NAME  THE  OUTPUT  FILE’) 

READ(*,10)OUNIT 

OPEN(8,FILE=OUNIT,MODE=’WRITE’ ,STATUS=’NEW’ ) 

GET  THE  INPUT 

C 

100  WRITE!*, 1) 

1  FORMAT!’  INPUT  THE  TARGET  LATITUDE,  LONGITUDE  AND  THE  RANGE’,/ 

1  ’  -  FREE  FIELD  -  RANGE  =  0  STOPS  EXECUTION’,/) 

READ!*,*)LAT,LON,LLR 

IF!LLR.LE.O)STOP 


non  ooo  non  noo  ooo  noo  ononoo 


STRETCH  OUT  IN  THE  LONGITUDINAL  DIRECTION 

LTR=LLR/COS ( FLOAT ( LAT ) /57 . 3 ) 

START  INPUT  UNIT  4  AT  THE  BEGINNING 

REWIND  4 
MS  =  0 

READ  A  HEADER 

200  READ ( 4 , 2 , END= 300 ) NST , I LAT , I  LON , I TM , NREC 

2  FORMAT (BZ,3X,I7,4X,2I5,8X,I2,I5) 

IF(ABS(ILAT-LAT).LE.LLR  .AND.  ABS( ILON-LON ) . LE. LTR )  THEN 

IT’S  IN  THE  RECTANGLE,  GET  THE  DATA 

MS=MS  +  1 
NSTAT(MS)=NST 
LATS(MS)=ILAT 
LONS(MS)=ILON 

READ(4,3)( (FEd ,MS,J) ,J=1,16),YRMD(I ,MS) ,I=1,NREC) 

3  F0RMAT(BZ,16F4.0,I8) 

ELSE 

OR  SKIP  IT 

READ(4,4)(JUNK, 1=1, NREC) 

ENDIF 

4  FORMAT(Il) 

GO  BACK  FOR  MORE  IF  WE  DON’T  HAVE  THE  MAX 

IF(MS.LT.NS)  GO  TO  200 
PRINT  7 

7  FORMAT!’  STOPPING  BEFORE  END’) 

300  CONTINUE 
IF(MS.LT.4)  THEN 

WRITE!*, 5)MS 

NOT  ENOUGH  DATA,  START  OVER 

5  FORMAT !'  ONLY ’,12,’  STATIONS  FOUND  -  START  OVER’) 

GO  TO  100 

ELSE 

WRITE!*,6)MS,!NSTAT!N),N=1,MS) 

6  FORMAT! 1 5,’  STATIONS  FOUND: ’ ,/( 1 10 ) ) 

ENDIF 

GET  THE  DATE  RANGE 

301  WRITE!*, 8) 


nno  ooo  non 


8  FORMATC  input  the  beginning  and  ending  days  0  =  1,  <  =  152  )V 

1  ’  BEGINNING  =  ZERO  GOES  TO  NEW  LATITUDE/LONGITUDE’/ 

2  ’  BEGINNING  =  NEGATIVE  ENDS  RUN’) 

READ(*,*)IDB,IDE 

IF( IDB.LT.O)STOP 
IF(IDB.EQ.O)GO  TO  100 
C 

C  OK,  LET’S  DO  THE  JOB 

C 

C  ZERO  OUT  FLAG  FIELDS 

C 

DO  304  L=l,16 
DO  303  N=1,NS 
NDLF(N,L)  =  0 
NDHF(N,L)  =  0 
NSDF(N,L)  =  -1 
NMNLF(N,L)  =  -1 
NMNHF(N,L)  =  -1 
NDG(N,L)  =  0 
NDDT(N,L)  =  0 

303  CONTINUE 

304  CONTINUE 

LOOP  OVER  ALL  LEVELS 

DO  400  L=l,16 
KS=0 

KS  COUNTS  STATIONS  WITH  A  HISTORY  OF  2  OR  MORE  DAYS 

DO  305  N=1,MS 
RNT(N)  =  0 
XMN(N)  =  0. 

KNT  COUNTS  HOW  MANY  DAYS 

305  CONTINUE 
DO  310  I=IDB,IDE 

ND=0 

DO  308  N=1,MS 

C 

C  ND  COUNTS  THE  NUMBER  OF  STATIONS  WITH  GOOD  DATA  TODAY 

C  FED  SAVES  IT  FOR  THE  TEST 

C  NSD  TELLS  WHICH  STATIONS  HAVE  GOOD  DATA 

C  FER  SAVES  IT  ALL  FOR  MEAN  AND  STD  CALCS 

C  YRMD  SAVES  THE  DATE  INFORMATION  (NOT  USED  PRESENTLY,  6/21 

C 

FED(N)=0. 

IFIFEd  ,N,L).LT.999)  THEN 
NDG(N,L)=NDG(N,L)+1 
KNT(N)  =  KNT(N)  +  1 
ND=ND+1 

FEU(ND)=FE( 1 ,N,L) 

NSD{ND)=N 


I  j. 


FER(KNT(N) ,N)  =  FE(I,N,L) 

YMD(KNT(N),N)  =  YRMD(I,N) 

XMN(N)=XMN(N)+FER(KNT(N),N) 

ENDIF 

308  CONTINUE 
C 

C  TEST  CURRENT  DAY’S  DATA 

C 

CALL  R11P05{FED,ND,NB,IB,NT,IT,IER) 

IF(IER.LE.O)  THEN 

IF{NB.GT.O)NDLF(NSD(IB),L)  =  NDLF(NSD(1B) ,L)  +  1 
IF(NT.GT.O)NDHF(NSD(IT),L)  =  NDHF( NSD( IT ) , L)  +  1 
DO  309  N=1,ND 

NDDT(NSD(N),L)  =  NDDT{NSD(N) ,L)  +  1 

309  CONTINUE 
ENDIF 

310  CONTINUE 

DO  315  N=1,MS 

IF(NDDT(N,L).EQ.O)  THEN 
NDLF(N,L)  =  -1 
NDHF(N,L)  =  -1 
ENDI F 

315  CONTINUE 

DO  330  N=1,MS 

IF(KNT(N).GT,  DTHEN 
C 

C  CALCULATE  MEAN  AND  ST.  DEV  IF  MORE  THAN  ONE  DAY  OF  DATA 

C 

KS=KS+1 

NSS(KS)  =  N 

XMN(KS)=XMN(N)/KNT(N) 

S=  0. 

DO  320  I=1,KNT(N) 

S=S+(FER(I ,N)-XMN(KS) )**2 
320  CONTINUE 

STD(KS)  =  SQRT(S/(KNT(N)-1)) 

ENDIF 

330  CONTINUE 
C 

C  RUN  TEST  ON  MEAN  VALUES  FOR  THE  STATIONS 

C 

CALL  R11P05(XMN,KS,NB,IB,NT,IT,IER) 

IF{IER.EQ.O)THEN 
DO  340  N=1,KS 

NMNLF{NSS(N) ,L)  =  0 
NMNHF{NSS(N) ,L)  =  0 
340  CONTINUE 

IFINB.NE.O)  NMNLF(NSS(IB),L)=1 
IF(NT.NE.O)  NMNHF(NSS(IT),L)=1 
ENDIF 
C 

C  RUN  TEST  ON  STANDARD  DEVIATIONS  OF  THE  STATIONS 

C 

CALL  R10P05(STD,KS,NT,IT,IER) 


on  o  o  o 


IF(IER.E0.0)THEN 
DO  360  N=1,KS 

NSDF{NSS{N),L)  =  0 
360  CONTINUE 

IF(NT.GT.O)NSDF{NSS(IT),L)  =  1 
ENDIF 

400  CONTINUE 
MDATB=0 
MDATE=0 
DO  410  N=1,MS 

I F ( MDATB . EQ . 0 ) MDATB= YRMD ( I DB , N ) 

I F ( MDATE . EQ . 0 ) MDATE=YRMD { I DE , N ) 

410  CONTINUE 

I F ( MDATB . EQ . 0 ) MDATB=MDATE- ( I DE- I DB ) 

I F { MDATB . LT . 0 ) MDATB= I DB 
I F ( MDATE . EQ . 0 )MDATE=MDATB+ ( I DE- I DB ) 

OUTPUT  RESULTS  FOR  VARIOUS  TESTS 

WRITE{8, 11 )MDATB, MDATE, ITM,LAT, LON, LLR.LTR 

11  FORMATCI  RESULTS  OF  THE  OUTLIER  TESTS  FROM’, 

1  18, ’  TO’ ,18, ’  AT’ ,13, ’Z’// 

2  ’  NUMBER  OF  FAILURES  BY  LEVEL  AND  STATION:  .=0,  *=NOTEST’ 

3  //’  TARGET  LATITUDE/LONGITUDE  IS’, 215,’,  RANGE: ’ , 1 2 , ’ X ’ , 1 2 

4  //’  RESULTS  OF  THE  DAILY  OUTLIER  TEST  (LOW  FAIL/HIGH  FAIL)’) 
WRITE(8,12)((L),L=1,16) 

12  FORMAT(/’  STATION  /  Lvl’,1614/) 

KRITE(8,13)(NSTAT(N), 

1  {SY(NDLF(N,L) ),SY(NDHF(N,L)),L=1,16),N=1,MS) 

13  FORMAT(I8,5X,32A2) 

WRITE(8,20)IDE-IDB+1 

20  FORMAT!//’  NUMBER  OF  DAYS  DAILY  OUTLIER  TEST  WAS  PERFORMED’ 

1  ’  OUT  OF’, 14) 

WRITE{8,12)( (L),L=1,16) 

WRITE(8,18)(NSTAT(N),(NDDT{N,L),L=1,16),N=1,MS) 

18  FORMAT! 18, 4X, 1614) 

WRITE(8,14) 

14  FORMAT!//’  RESULTS  OF  THE  MEANS  OUTLIER  TEST  (LOW  FAIL/HIGH’ 

1  ’FAIL)’) 

WRITE!8,12)!!L),L=1,16) 

WRITE!8,13)!NSTAT!N)  , 

1  !SY!NMNLF!N,L)),SY!NMNHF(N,L)),L=1,16),N=1,MS) 

WRITE(8,15) 

15  FORMAT!//’  RESULTS  OF  THE  STANDARD  DEVIATION  OUTLIER  TEST’) 
WRITE!8,12)!!L),L=1,16) 

WRITE(8,16)!NSTAT!N),!SY!NSDF!N,L)),L=1,16),N=1,MS) 

16  FORMAT!I8,4X,16A4) 

WRITE!8,17)IDE-IDB+1 

17  FORMAT!//’  NUMBER  OF  DAYS  REPORTING  OUT  OF ’,14) 
WRITE!8,12)!!L),L=1,16) 

WRITE!8,I8)!NSTAT!N) ,!NDG!N,L),L=l,16),N=l,MS) 

OUTPUT  SUMMARY  BY  STATION 


i  b 


c 

c 

c 

c 

c 

c 

c 


WRITE(8,22) 

22  FORMAT (//’  ************’, 

1  't**t1lit**********************************^^******^****1^1^t***t.*1^**  * 

2  //’  SUMMARY  BY  STATION  FOLLOWS’/) 

DO  500  N=1,MS 

WRITE! 8,23 )NSTAT(N ) ,LATS(N),LONS(N) 

23  FORMAT!//’  - ’, 

1' - 

2  /’  STATION’ ,16, ’ :  LAT.LON  =  ’,14, ’,’,14/ 

3  /’  LEVEL  #DAYS  DAILY  TESTS  MEANS  TESTS 

4  ’  STD  TESTS  TOTAL’/) 

DO  420  L=l,16 

SET  CODES  FOR  TESTS:  -1  IS  NOTEST,  0  IS  PASS  ALL,  1  IS  FAIL 

ONE  OR  MORE 

FOR  THE  TOTAL  TEST,  NOTEST  MEANS  NO  TESTS  PERFORMED, 

PASS  MEANS  PASS  ALL  TESTS  PERFORMED  !IF  ANY) 
FAIL  MEANS  FAILED  AT  LEAST  ONE 

NFL=MIN!NDLF!N,L) ,  1 ) 

NFH=M1N!NDHF!N,L),1) 

NMFL=MIN!NMNLF!N,L),1) 

NMFH=MIN!NMNHF!N,L),1) 

NSF=MIN!NSDF!N,L),1) 

NTF=MAX ! NFL , NFH , NMFL , NMFH , NSF ) 

WR1TE!8,24)L,NDG!N,L) ,FLAGSL!NFL),FLAGSH!NFH), 

1  FLAGSL ! NMFL ) , FLAGSH ! NMFH ) , FLAGS ( NSF ) , FLAGS ( NTF ) 

420  CONTINUE 

24  FORMAT !  2 1 6 , 4X ,  A6 , 3X ,  A6 , 5X ,  A6 , 3X ,  A6 , 5X ,  A6 , 5X ,  A6  ) 

500  CONTINUE 

WRITE!8, 19) 

19  FORMAT!//'  ******************’ ^ 

1  ’  iit**t:*******t*t******************************************  ’/////) 

GO  TO  301 
END 


SUBROUTINE  R10P05 ! X , N , NOUTT , IT, I ER ) 

PARAMETER  !MAX=30) 

DIMENSION  X!N),SX!MAX),ISX!MAX) 

DIMENSION  CV!30) 

C 

C  THIS  SUBROUTINE  DOES  THE  RIO  TEST  FOR  OUTLIERS  AT  THE  .05  CRITICAL 
C  TEST  LEVEL. 

C 

C  THE  CV  VALUES  ARE  THE  CRITICAL  VALUES  FROM  TABLE  V. 

C 

DATA  CV/2*0, .941, .765, .642,. 560,. 507,. 468, .437, .412, .392, .376, 

1  . 361 , . 349 , . 338 , . 329 , . 320 , . 313 , . 306 , . 300 , . 295 , . 290 , . 285 , . 281 , 

2  .277, .273, .269, .266, .263, .260/ 

C 

C  INPUT  VALUES  ARE:  X  -  THE  DATA  TO  BE  TESTED 
C  N  -  NUMBER  OF  DATA  POINTS 

C 

C  OUTPUT  VALUES  ARE:  NOUTT  -  FLAG  FOR  LARGEST  VALUE 


oononnoon  oonnn 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 


0  IF  OK,  1  IF  NOT 
IB  -  LOCATION  IN  X  ARRAY 

lER  -  ERROR  RETURN  INDICATOR 

0  IF  ALL  OK 

1  IF  N  OUT  OF  RANGE  4  TO  MAX,  INCLUSIV 
NO  TEST  PERFORMED 


CHECK  FOR  SUITABLE  N 


lER  =  0 
NOUTT  =  0 

IF(N.LT.3.0R.N.GT.MAX)THEN 
lER  =  1 
RETURN 
END  IF 

DO  100  I=1,N 
SX(I)  =  Xd) 

ISXd  )  =  I 
100  CONTINUE 


C 

C  SORT  INTO  INCREASING  ORDER 
C 

CALL  SHSORT(SXdSX,N) 

RT  =  0. 

IF(SX(N)-SX{1)  ,NE.  0.)  RT  =  (SX(N)  -  SX(N-1 ) )/(SX(N)  -  SX(1)) 
IF{RT.GT.CV'(N)  )THEN 
NOUTT  =  1 
IT  =  ISX(N) 

ENUl  F 

RETURN 

END 


SUBROUTINE  R I 1P05 { X , N , NOUTB , I B, NOUTT , I T , I ER ) 

DIMENSION  X(N),SX(30),ISX(30) 

DIMENSION  CV(30) 

THIS  SUBROUTINE  DOES  THE  Rll  TEST  FOR  OUTLIERS  AT  THE  .05  CRITICAL 
TEST  LEVEL. 


THE  CV  VALUES  ARE  THE  CRITICAL  VALUES  FROM  TABLE  V. 

DATA  CV73*0, .955, .807, .689,. 610,. 554,. 512, .477, .450, .428, .410, 
1  .395, .381 , .369, .359, .349, .341,. 334, .327,. 320, .314, .309, .304, 

1  .299, .295, .291, .287, .283/ 


INPUT  VALUES  ARE:  X  -  THE  DATA  TO  BE  TESTED 

N  -  NUMBER  OF  DATA  POINTS 

OUTPUT  VALUES  ARE:  NOUTB  -  FLAG  FOR  SMALLEST  VALUE 

0  IF  OK,  1  IF  NOT 
IB  -  LOCATION  IN  X  ARRAY 
NOUTT  -  FLAG  FOR  LARGEST  VALUE 
0  IF  OK,  1  IF  NOT 


1  7 


o  o  o  n  o 


C 

C 

c 

c 

c 

c 

c 

c 


c 

c 

c 


IT  -  LOCATION  IN  X  ARRAY 

lER  -  ERROR  RETURN  INDICATOR 

0  IF  ALL  OK 

1  IF  N  OUT  OF  RANGE  5  TO  30,  INCLUSIVE 


CHECK  FOR  SUITABLE  N 
lER  =  0 

IF(N.LT.4.0R.N.GT.30)THEN 
lER  =  1 
RETURN 
ENDIF 

DO  100  1=1, N 
SX(l)  =  X(I) 

ISXd  )  =  I 
100  CONTINUE 

SORT  INTO  INCREASING  ORDER 

CALL  SHSORT(SX,ISX,N) 

RB  =  0. 

IF(SX(N-1 )-SX{ 1 )  .NE.  0.)  RB  =  (SX(2)  -  SX(  1 ) )/(SX(N-l )  -  SX(1)) 
NOUTB  =  0 

IF(RB.GT.CV(N) )THEN 
NOUTB  =  1 
IB  =  ISX(l) 

ENDIF 
RT  =  0. 

IF(SX(N)-SX(2)  .NE.  0.)  RT  =  (SX(N)  -  SXIN-I) )/(SX(N)  -  SX(2)) 
NOUTT  =  0 

IF(RT.GT.CV(N) )THEN 
NOUTT  =  1 
IT  =  ISX(N) 

ENDIF 

RETURN 


END 


C  Included  for  completeness  only  -  from  NONIMSL  library  at  NFS 


C  . 

C 

C  A.  idem  I  fra  11  ON: 
r  TITLE:  NUMERICAL  SORT 

C  ID:  Ml-NPG-SHSORT  (F-IV) 

C  PROGRAMMER:  R,  BRUNELL 

C  DATE:  MARCH  1968 

C  MODIFIED:  DEC.  1973  BY  L.  NOLAN 

B.  PURPOSE: 

TO  SORT,  IN  ASCENDING  ORDER,  AN  ARRAY  OF  SINGLE  PRECISION  REAL 
NUMBERS  BY  THE  METHOD  OF  SHELL,  AND  TO  PRODUCE  AN  ARRAY  OF  INDEXES 
SO  USER  CAN  RE-ORDER  OTHER  CORRESPONDING  INFORMATION  ACCORDING  TO 
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ASCENDING  VALUES  OF  "A". 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 


C.  USAGE: 

1.  CALLING  STATEMENT: 

CALL  SHSORT(A,KEY,N) 

2.  ARGUMENTS: 

A  -  ARRAY  OF  NUMBERS  TO  BE  SORTED.  THIS  ARRAY  IS  SORTED 
(RE-ORDERED)  BY  "SHSORT". 

KEY  -  ARRAY,  DIMENSIONED  AT  LEAST  N  IN  CALLING  PROGRAM,  TO  BE 
FILLED  BY  USER  WITH  INTEGERS  FROM  1  TO  N.  AFTER  EXIT 
FROM  SHSORT,  KEY(l)  WILL  CONTAIN  THE  ORIGINAL  INDEX  OF 
THE  SMALLEST  ELEMENT  OF  "A";  KEY(2)  WILL  CONTAIN  THE 
ORIGINAL  INDEX  OF  THE  NEXT-TO-SMALLEST  ELEMENT  OF  "A"; 
ETC.  KEY(N)  WILL  CONTAIN  THE  ORIGINAL  INDEX  OF  THE 
LARGEST  ELEMENT  OF  "A”. 

N  -  NUMBER  OF  MEMBERS  IN  ARRAYS  "A"  AND  "KEY". 

D.  REFERENCES: 

1.  "ALGORITHM  201,  SHELLSORT" ,  BOOTHROYD,  J.,  "COMMUNICATIONS  OF 

ACM",  VOL  6,  NO  8,  AUGUST  1963,  P.445. 

2.  "CERTIFICATION  OF  ALGORITHM  201",  BATTY, M. A..  "COMMUNICATIONS 

OF  ACM",  VOL  7,  NO  6,  JUNE  1964,  P.349. 


SUBROUTINE  SHSORT (A, KEY, N) 
DIMENSION  A(N),KEY(N) 

Ml  =  l 

6  M1=M1*2 

IF  (Mi  .LE.  N)  GO  TO  6 

Ml=Ml/2-l 

MM=MAX0(M1/2,1) 

GO  TO  21 

20  MM=MM/2 

IF  (MM  .LE.  0)  GO  TO  100 

21  K=N-MM 

22  DO  1  J=1,K 
II=J 

11  IM=II+MM 

IF  (A(IM)  .GE.  A(II))  GO  TO  1 
:emp=a( II  ) 

IT=KEY(II ) 

A(II  )=A(IM) 

KEYdl  )=KEY(IM) 

A(IM)=TEMP 
KEY( IM)=IT 
1 1=1 1 -MM 

IF  (II  .GT.  0)  GO  TO  11 
1  CONTINUE 
GO  TO  20 
100  RETURN 
END 
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